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Abstract 

We calculate the thermomechanical properties of a-iron, and in particular its isothermal and 
adiabatic elastic constants, using first-principles total-energy and lattice-dynamics calculations, 
minimizing the quasi-harmonic vibrational free energy under finite strain deformations. Particular 
care is made in the fitting procedure for the static and temperature-dependent contributions to 
the free energy, in discussing error propagation for the two contributions separately, and in the 
verification and validation of pseudopotential and all-electron calculations. We find that the zero- 
temperature mechanical properties are sensitive to the details of the calculation strategy employed, 
and common semi-local exchange-correlation functionals provide only fair to good agreement with 
experimental elastic constants, while their temperature dependence is in excellent agreement with 
experiments in a wide range of temperature almost up to the Curie transition. 
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I. INTRODUCTION 


Elemental iron is a material of great scientific and economic interest: it’s the major 
constituent of steels, it determines the properties of the earth core, and its complex phase 
diagram is driven by the subtle interplay between vibrational and magnetic contributions, 
making it particularly challenging to describe accurately from hrst-principles. This is espe¬ 
cially relevant as temperature increases, since magnetic excitations become important and 
a dramatic change in the magnetic nature of the system takes place. At ordinary pressures, 
iron turns from a BCC ferromagnet to a BCC paramagnet with a second-order transition 
{a —)■ /3) at a Curie temperature of ~ 1043 K. This transition is then followed by two struc¬ 
tural transitions, a BCC—?-FCC (/3 —)■ 7 ) transition at 1185 K and a FCC—?-BCC (7 —)■ 6) 
transition at 1670 K, before hnally melting at ~ 1814 K. 

First-principles simulations can be a key technique and a valid alternative to experiments 
in order to get accurate predictions of phase diagrams without the need of phenomeno¬ 
logical parameters, and become essential at conditions that are challenging to reproduce 
in real life, like those inside the earth’s cor^. For the case of pressure-temperature phase 
diagrams, zero-temperature hrst-principles equations of state can be supplemented with 
hnite-temperature vibrational entropies, that can be derived directly from the knowledge of 
the phonon dispersions. These latter can be calculated from hnite differences, or more ele¬ 
gantly and less expensively with density-functional perturbation theory (DFPT)I^. When 
coupled to the quasi-harmonic approximationP (QHA), these techniques allow to calculate 
thermal expansion and vibrational properties at hnite temperatures, often well above the De¬ 
bye temperaturd^^. While there have been numerous hrst-principles calculations of elastic 
properties of solids by either total energy, stress-straii J^^hsi qj. density-functional perturba¬ 
tion theory approached®® a relatively limited number of them has been focused on to the 
thermo-mechanical properties of metals or minerald®®. 

In this paper, we calculate the adiabatic and isothermal hnite-temperature elastic proper¬ 
ties of ferromagnetic a-iron at ambient pressure and in the temperature range of stability for 
ferromagnetic a-iron using the QHA and DFPT. We also carefully explore multidimensional 
htting procedures for the static and vibrational contributions to the free energy, and analyze 
the quality of the ht and the source of errors of both contributions, providing a conhdence 
interval of each elastic constant as a function of temperature. 
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The paper is organized as follows: in Sec. |TT] we introduce the finite strain framework 


used to calculate the elastic constants, and in Sec. we give the computational details 
of our first-principles density-functional theory (DFT) and DFTP calculations. We present 


our results and their comparison to experiments in Sec. [IVj Finally, Section |V] is devoted to 
summary and conclusions. 


II. FINITE-STRAIN METHOD 

In the limit of small deformations, the energy of a crystal at an arbitrary conhguration can 
be Taylor-expanded in terms of a symmetric matrix e describing a uniform linear deformation 
A such that 


A = 1 + e 


( 1 ) 


and any position vector f in the reference conhguration is transformed into (1 -|- s) ■ f. The 
isothermal elastic constants (SOECs) at zero pressure are then dehned as the second-order 
coefficients of this expansion according tcP^ 
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where F is the Helmholtz free-energy and £*, Sj,i,j = 1...6 are the components of the 
strain tensor £ (we adopt here the Voigt notation). Also, note that the second derivatives are 
evaluated at the thermodynamic equilibrium conhguration with volume Vq and at constant 
temperature T. 

For cubic crystals, like a-iron, only three elastic constants are needed to completely 
determine the stihness tensor and, therefore, fully characterize the mechanical response of the 
system in the linear elastic regime. As a consequence only three independent deformations 
are sufficient, and we choose here the hydrostatic, tetragonal and trigonal deformations, 
shown in Tab. Each deformation fully determines one of the cubic elastic constants (or 
elastic moduli). 

In order to compute hnite-temperature properties and, therefore, to calculate the 
Helmholtz free energy F, the vibrational contributions must be added to the static energy 
contributions. The QHA® provides an analytical expression for the vibrational contribution 
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TABLE I. Deformations and corresponding strain vectors in the Voigt notation: hydrostatic, 

^^Hetragonal and ^^Hrigonal deformations are governed by a single parameter. *Note that the 
trigonal deformation reported here is the first-order expansion of the full strain tensor with 
non-zero off-diagonal terms. 

to the free energy: 
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where the sum is performed over all the phonon modes A and all the phonon wave vectors 
q spanning the Brillouin zone (BZ). Here, ks is the Boltzmann constant and Uq^x the vibra¬ 
tional frequencies of the different phonon modes, where in the QHA their explicit dependence 
on the geometry of the system via the primitive lattice vectors {oj} is accounted for. The 
vibrational part, coming directly from the analytic partition function of a Bose-Einstein gas 
of harmonic oscillators, can be split into a zero-point energy term plus a contribution which 
depends explicitly on the temperature T. We neglect here the thermal electronic effects, 
because they are believed to be small compared to the quasi-harmonic vibrational contribu- 
tiorpIEII in range of stability of the a phase. Magnetic effects are also not considered, 
except for the longitudinal relaxation of the total magnetic moment as a function of strain, 
but they are known to be important approaching the Curie poinlP^^^ and their influence on 


the elastic properties is briefly discussed in Sec. |IV C| in the light of previous studies. 

Eq. is used to calculate the isothermal elastic constants at hnite temperature; however, 
in order to compare results with experimental data obtained from ultrasonic measurement^. 
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we also calculate the adiabatic elastic constants, using the following relation^ 
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where the heat capacity at constant volume Gy and the volumetric thermal expansion coeffi¬ 
cient a are both calculated from the QHA. The superscripts (S) and (T) dehne the adiabatic 
and isothermal elastic constants Gij and bulk modulus B respectively. 


III. COMPUTATIONAL DETAILS AND PSEUDOPOTENTIAL SELECTION 

We calculate the hrst-principles elastic constants using DFT, as implemented in the PWSCF 
and PHONON packages of the Quantum-ESPRESSO distribution!23 for the static and lattice 
dynamical calculations, respectively. The calculations are spin-polarized and the magnetic 
moment is free to vary collinearly in order to minimize the total energy. In all calculations the 
exchange-correlation effects have been treated within the generalized-gradient approximation 
(GGA) with the PBE functionaP^. We use an ultrasoft pseudopotentiaP^ (USPP) from 
pslibrary.0.3.(P^, which includes also 3s and 3p semicore stated (i.e. 16 valence electrons) 
along with a plane-wave basis with a wavefunction kinetic-energy cutoff of 90 Ry and a 
cutoff of 1080 Ry for the charge density. We sampled the BZ with an offset 24 x 24 x 24 
Monkhorst-Pack fc-mesh, with a Marzari-Vanderbilt smearin^^ of 0.005 Ry. 

Phonon calculations were carried out for each deformation within DFPT® the dynamical 
matrix and its eigenvalues are calculated on a 4 x 4 x 4 mesh of special points in the BZ and 
Fourier-interpolated on an extended 21 x 21 x 21 grid for the integration of thermodynamic 
quantities. We arrived at this computational setup (cutoff, smearing and BZ sampling) after 
a careful investigation of the convergence of total energy and individual phonon frequencies 
for different deformations. Also, we verihed that individual total energies and phonon fre¬ 
quencies do change smoothly as a function of strain. 

Since the choice of the pseudopotential is of primary importance for a clear comparison 
with computational and experimental data in the literature, it is worth to stress that the 
one used here has been chosen among different candidates from the pslibrarj^ and GBRV 
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FIG. 1. Equilibrium lattice parameter at 0 K for the different iron pseudopotentials tested in this 
work. All the data shown here are obtained with the PBE XC functional except for the last three 
columns on the right, where we have used PBE+tP^Hlol^ WCp^and a PBEsoP^ respectively (here 
we use a Hubbard U correction with U = 3eV). The data come from a Birch-Murnaghan fit, 
do not include zero-point energy and are compared to all-electron WIEN2Kp31, exciting® and 
VASF® calculations from Ref. I451I17I respectively and experiment d^^ * ^^ *^ (horizontal yellow line). 
The crosshatch dotted column corresponds to the pseudopotential chosen for the production runs. 


library® to reproduce, as closely as possible, the all-electron FLAPW equilibrium lattice 
parameter, bulk modulus at 0 K and local magnetization obtained from independent groups. 
Also, for the sake of completeness, we compare against results obtained using the VASP 
code and associated pseudopotentiald^. 

These values in Figs. 01 are obtained from a Birch-Murnaghan £t of calculated E(V) 
data points. Interestingly, we have found that the volume range of validity for htting a Birch- 
Murnaghan curve is limited on the expansion side due to anomalies in the E{V) curve and 
its derivatives. These anomalies, also reported for all-electron and other calculation methods 
in Ref. [52], are more clearly visible as “shoulders” in the M(V) behavior (see Fig. and, as 
visible from Fig. can be associated to a smooth magnetic transition from a low to high 
spin state due to the splitting of the majority and minority spin t 2 g electrons upon increasing 
the volume. However, for the pseudopotential chosen here, the expanded volumes at which 
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FIG. 2. Equilibrium bulk moduli at 0 K for the different iron pseudopotentials tested in this 
work. All the data shown here are obtained with the PBE XC functional except for the last three 
columns on the right, where we have used PBE+tPZHini^ WCP^ and a PBEsoP^ respectively (we 
use here a Hubbard U correction with U = 3eV). The data come from a Birch-Murnaghan fit, 
do not include zero-point energy and are compared to all-electron WIEN2Kp31, excitin^^ and 
VAS recalculations from Ref. USHSl respectively and experiment^e (horizontal yellow line). The 
crosshatch dotted column corresponds to the pseudopotential chosen for the production runs. 


this anomaly is observed (above 99{!^ are far beyond the theoretical thermal expansion of 
the system in the thermodynamic region considered in this work, thus enabling us to £t the 
energy surface with volume expansions up to ~9% still using a standard Birch-Murnaghan 
equation. 
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FIG. 3. (Top panel) Equation of state as a function of percent volume change with respect to the 
theoretical equilibrium configurations for three of the selected PBE pseudopotentials considered 
in this work (circles). The yellow circles best match the all-electron WIEN2k®I (pentagons) and 
excitin^SI (triangles) results from Ref. 0^1 and ED and correspond to the rrkjus-0.2. l-16e 
pseudopotential used in this work. Continuous lines are the best fit of the Birch-Murnaghan 
equation. (Bottom panel) Total magnetization as a function of percent volume change. The soft 
magnetic transition discussed in the text is visible as a clear change in the average slope of the 


different curves. 
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FIG. 4. [a] DOS of majority/minority (red/blue) spin channels at the equilibrium (solid line) 

and AV ~11% (dashed line) where, for the pseudopotential used in production run, the magnetic 
transition takes place, [b] The contribution of the t 2 g electrons to the majority/minority DOS 
(green/black) is also reported. To obtain a smooth DOS, a non-self-consistent calculation with an 
offset 60 X 60 X 60 Monkhorst-Pack A:-mesh is performed on top of a scf loop. 
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IV. RESULTS 


In this section we present results for selected thermodynamic quantities and for the three 
strain deformations (hydrostatic or volumetric, tetragonal, and trigonal). Each deformation 
determines uniquely one of the three elastic constants: B (bulk modulus), Cn and U 44 
respectively. 


A. B — Volumetric strain 

The volumetric deformation ed) pg described by a single parameter Sa, namely, the 
strain of the cubic lattice parameter (see Tab. |^. Thus, the lattice spacing is dehned as: 

a = ao(l + Ea) (5) 

where oq is the theoretical equilibrium lattice parameter without zero-point contribution (see 
Tab. III). The static part of the Helmholtz free energy of Sec.j^is obtained by htting a Birch- 
Murnaghan equation of statd^ to a series of well converged total energy values calculated 
on a one dimensional regular grid with Sa going from —0.02 to -1-0.03 in steps of 0.001. The 


resulting static contribution to the bulk modulus is reported in Tab. Ill The vibrational 


contribution, on the other hand, has been calculated on a coarser grid via integration of 
the phonon dispersions as from Eq. (examples for the calculated phonon dispersion and 
resulting Griineisen parameters can be found in Fig. [^, with Ea ranging from —0.012 to 
-1-0.020 in increments of 0.004 and htted with a second-order polynomial as a function of the 
strain parameter Ea- The stability of the results has been checked against a £t with lower 
and higher order polynomials (see Supplemental Material^). 

The free energy is then obtained as an analytical function of Ea and T and is shown in 
Fig. 1^ We then determined the thermal expansion (Fig. |^, the thermal expansion coefh- 
cients (Fig. [^, the heat capacity (Fig. and the isothermal bulk modulus from the 

analytic second derivative of the free energy as in Eq. The adiabatic correction of Eq. 
is then used to compute the adiabatic bulk modulus B^^\T). Results are reported in Fig. 
and compared to experimental data from Refs. [26] and [501 The agreement between experi¬ 
ments and calculations in the thermal behavior of the bulk modulus is remarkable, especially 
below the Debye temperature (0/) — 500 K). Above ©d, the small deviation from experi¬ 
ments can be ascribed to magnetic fluctuation^^SEMSEll become increasingly important 
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FIG. 5. (Left panel) Phonon dispersions along high-symmetry directions in the BZ calculated at the 
theoretical electronic equilibrium volume (blue solid line) and at the quasi-harmonic theoretical 
equilibrium volume at 300 K (red dashed line). The results (see also Fig. 1 in Supplemental 
Material and Ref. |8] or Ref. [56] for comparison with previous theoretical data) are compared to 
experimental data at room temperature from Ref. [57| (Expt.l - squares) and Ref. [58] (Expt.2 - 
circles). (Right panel) Griineisen parameters calculated along the same path in the BZ and the 
same equilibrium volumes used for the phonon dispersion (blue solid line for the 0 K case and 
red dashed line for the 300 K case). The Griineisen parameters are obtained computing the first 
derivative with respect to the volume of a cubic fit of the phonon frequencies. 

See EPAPS Document No.[] 


approaching the Curie temperature (1043 K), plus minor contributions from anharmonic 
effects (beyond quasi-harmonic) and from the electronic entropy. At 1000 K, the softening 
of the calculated is nearly 15% with respect 0 K. The calculated magnetic moment in¬ 
creases from 2.17 /i^ per atom (2.22 from experiments^ at the 0 K equilibrium volume 
to 2.27 at the 1000 K equilibrium volume. Obviously, transverse magnetic fluctuations 


are neglected in these calculations, and we postpone to Sec. jlV Cj the discussion on the 
mismatch between experiments and calculations in absolute values. 


B. Cii, (744 ~ Tetragonal and Trigonal strains 

The Helmholtz free energy F depends upon two strain parameters: the isotropic lattice 
strain Sa and a second strain parameter or Ed according to the deformation considered 
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FIG. 6. (Top panel) Free-energy landscape of cubic BCC iron as a function of volume V and 
temperature T. The dashed black line corresponds to the set of points that minimize the free- 
energy surface at each temperature. The continuous green and blue lines are the projections of the 
black dashed line in the T-V and F-T planes, thus describing the volumetric thermal expansion 
and the zero-pressure free energy as a function of T. (Bottom panel) Volumetric thermal expansion 
(green solid line) compared to experimental data from Ref. [36] (Expt.l ~ blue squares, note that 
below room temperature the data are extrapolated according to the thermal expansion coefficient 
of Ref. and Ref. [59] (Expt.2 - magenta triangles). As a guide to the eye, we also report in all 
plots shifted and scaled quantities. The former are rigidly translated on the vertical axis while the 
latter are multiplied by a constant factor to match the experimental 0 K value. 

(see Tab. [^. 

The tensor is associated to a continuous tetragonal deformation that stretches the 
edge c of the cubic undistorted structure along the 2 ; axis while leaving unchanged the other 


12 












Temperature (K) 


FIG. 7. (Left panel) Linear (green solid line) coefficient of thermal expansion compared to ex¬ 
perimental data from Ref. [60] (squares). (Right panel) Specific heat at constant pressure (green 
solid line), at constant volume (blue dashed line), and compared to experimental data from Ref. 1611 
(Expt.l ~ squares) and from Ref. [62| (Expt.2 - triangles) . 



EIG. 8. Adiabatic bulk modulus as a function of T (blue continuous line) calculated along with 
its confidence interval on the fit (shaded green) and compared to experimental data from Ref. (501 
(Expt.l - green circles) and from Ref. |26| (Expt.2 - yellow triangles). As a guide to the eye, we 
also plot the bulk modulus rigidly shifted (dotted line) and scaled (dashed line) to match the 
experimental 0 K value. 

edges. The relation between the strain Sc and the distorted edge c is: 

c = a{l + ec). (6) 

The tensor is associated to a continuous trigonal deformation that stretches the main 
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diagonal d of the undistorted cubic structure along the (111) direction while tilting the 
undistorted edges and preserving their length. In this case, the relation between the strain 
Ed and the distorted main diagonal is 

d=V3a{l + ed), (7) 


while the relation with the cosine of the angle between the distorted edges is 


cos {a) 


1 ~ ^d(2 + Ed) 

i^d — l)(^d + 3) 


( 8 ) 


Both deformations do not conserve the volume per atom. In particular, in the tetragonal 
one the volume increases as a function of Ec, while in the trigonal case, the volume decreases 
as a function of Ed- Alternatively, we could have chosen volume-conserving deformations as 
in Ref. dH but the advantage of the present scheme is that each deformation determines 
uniquely one elastic constant at the time, and enables us to determine easily the conhdence 
interval of each elastic constant by error-propagation theory. 

In the next sub-sections we describe the calculation of the static and vibrational contri¬ 
butions, separately. The reason is that we want to analyze their contributions to the global 
energy landscape separately. This also allows us to sample the two contribution landscapes 
with two different grids. Indeed, the static term displays a minimum as a function of the 
strain parameters and has to be sampled with a dense grid, while, on the other hand, the 
vibrational term is flat, monotonic and can be sampled with a coarse grid. 


1. Static contribution 

To evaluate the static contribution to the elastic constants, we performed a series of well 
converged total energy calculations on a two dimensional discrete grid [Ea-,Sc/d\ (see Fig. 
for details on the grid). The Ea grid is asymmetric with respect to zero and with more points 
in the positive range of the strain parameter, in order to sample accurately the values of the 
static contribution to the free energy also in the thermal expansion range. 

The resulting total energies are htted with a two-dimensional bivariate polynomial up to 
5th degree using a least-square methocP^. 

The analysis of the quality of the £t of discrete data points to a two-dimensional energy 
surface is crucial to resolve the possible sources of error that could affect our elastic constants 
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Order AAE (Ry) 

tetragonal 

AAE (Ry) 

trigonal 

2 

1.894 10"® 0.997259 

9.163 10-5 0.971580 

3 

1.997 10-® 0.999975 

7.369 10-5 0.999819 

4 

1.002 10-® 0.999993 

2.933 10-5 0.999968 

5 

9.150 10-'^ 0.999995 

1.693 10-5 0.999990 


TABLE II. Average absolute error (AAE) and adjusted coefficient of determination (i?^) of the 
two-dimensional fit of the static energy landscape, for the tetragonal and trigonal deformations, as 
a function of the order of the polynomial. 


and, therefore, for a reliable comparison with experiments and the wide range of scattered 
data available in the literatnre. Therefore, in addition to the visnal inspection of the £t 
along constant sections, we evalnated the adjusted coefficient of determination (i?^) and 
the the average absolute error (AAE), defined as: 


AAE ^ ^ 




N ^ 


(9) 


where N is the total number of [ea,Sc/d\ discrete values and is the bivariate polynomial 
of degree n. Thus, is a measure of the quality of the fitting model, i.e. how well 
the analytic function approximates the calculated data points. The AAE is a quantitative 
measure of the distance of the fitted curve from the calculated points. We found that the 
AAE decreases by increasing the degree n of the polynomial and R? approaches unity, as 
shown in Tab. [TT| According to these results, in both cases, we considered the 4th-degree 
polynomial to provide a sufficiently accurate fit (indeed the AAE is two orders of magnitude 
smaller than the difference between the highest and the lowest total-energy data points). 

Fig. shows a plot of the static energy landscape for both the tetragonal and trigonal 
deformations, with the minimum elongated along the diagonal in the [sa, £c] space or along 
constant Sd in the space. 


2. Vibrational contribution 

In order to evaluate the vibrational contributions to the free energy, we performed a series 
of linear-response phonon calculations on a two dimensional grid in the space of deformation 
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FIG. 9. Static energy landscape of the tetragonal (top panel) and trigonal (bottom panel) distor¬ 
tions projected on the [ea-,£c/d\ plane. 


scalars Sa, Sc/d- Since the lattice dynamics calculations are one order of magnitude more 


time consuming than the total energy calculations, we used a coarser grid (see Fig. 10 for 
details on the grid). 

The eigenvalues of each dynamical matrix are Fourier-interpolated in order to obtain 
smooth and continuous phonon dispersions. The zero-point energy and the thermal contri¬ 
butions are calculated by numerical integration over 21x21x21 points in reciprocal space. 
This is essential to obtain numerically accurate values of the vibrational contribution. 

Like for the case of the static contribution, we determined the best polynomial necessary 
to £t our data over the entire temperature range from 0 to 1000 K. Similarly, we used the 
adjusted and the AAE as indicators of the quality of the £t. We also checked a posteriori 
the convergence of the elastic constant curves obtained by htting to different polynomial 
degrees. In the tetragonal case, a quadratic bivariate polynomial (i.e. 6 parameters) is 
sufficient to accurately reproduce the distribution of data points. On the other hand, for 
the trigonal deformation, a 4th order bivariate polynomial (i.e. 16 parameters) is needed. 
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FIG. 10. Vibrational quasi-harmonic contribution to the Helmholtz free energy at T = 750 K 
in the [eajGc] space (top panel), [£a,£d\ space (bottom panel). A 2nd and a 4th order bivariate 
polynomial are respectively used to fit the tetragonal and trigonal data sets. 


Our choice of polynomial is dictated by the need to minimize the AAE, maximize and 
minimize the conhdence interval as a function of temperature (see Supplemental MateriaP^ 
for the stability of the results against other polynomials). As an illustration, we report the 
vibrational energy landscape at 750 K for the tetragonal and for the trigonal distortions 


(Fig. 10). 
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FIG. 11. Free energy landscape of the tetragonal (top panel) and trigonal (bottom panel) 
tions projected on the and [£a-,£d\ plane respectively at 500 K. The green squares 




3. Evaluation of the elastic constants 


Next, we sum the static and vibrational energy landscapes obtained in 
sections and compute the Helmholtz free energy. An example of the resulting 


500 K is displayed in Fig. IT The second derivative with respect to strain can 
analytically at the minimum of the free energy as a function of temperature. 


the previous 
landscape at 
be evaluated 


Then, in order to understand if the discrepancy between the experimental and calculated 
elastic constants could be ascribed to the fitting procedure, we have calculated the confidence 
interval of Cn and 0^4^. To this end, we have computed the covariance matrix of each best-fit 
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contribution to the free energy, defined as: 


cov(P) = cr^ {J^J) ^) 


( 10 ) 


where P is the set of polynomial coefficients, is the squared residual and J is the Jacobian 
matrix which is provided in output by the least squares routine. The global variance of each 
best fit polynomial is then obtained by considering both the diagonal and the off-diagonal 
elements of the covariance matrix cov(P). Finally, we used error-propagation theory to 
obtain the confidence interval of the elastic constants. 


The calculated Cu and 6*44 elastic constants of BCC a-iron both decrease by increasing 


temperature, as shown in Fig. 1^ Our results are in reasonable accordance with those 
reported in Ref. [ 68 ] (the exception is 044 that in our case is fairly underestimated) where, 
however, a direct detailed comparison with experimental thermal softening is clearly more 


difficult. In Tab. Ill, we report their zero temperature values with and without zero-point 


energy (ZPE) contributions, thus comparing these to experimental values. Also, for sake 

i^the 012 = C = \{Cii - C 


of completeness, we report in Tab. 


12 , 


= !(C^: 


11 


B) 


and anisotropy ratio 04 , 4 ,/C obtained from standard theory of elasticity and our calculated 


P, Cu and C 44 (see Fig. 13 for the temperature dependence of the C). The inclusion of 
ZPE results in a small decrease of the elastic constants and bulk modulus. The confidence 
interval at zero temperature is of the order of 0.1 GPa and cannot account for the difference 
with respect to experiments, so we will discuss other possible source of this discrepancy in 
the next section. 


r(K) 

a (A) B (GPa) 

Cu (GPa) C 44 (GPa) 

0 (no ZPE) 

2.834 

199.8T0.1 

296.7T0.3 

104.7T0.1 

0 (ZPE) 

2.839 

194.6T0.3 287.9T0.4 

102.2T0.5 

0 (Expt.)P5ni 2.856 

170.3T1 

239.5±1 

120.7T0.1 


TABLE III. Calculated 0 K elastic constants for iron with and without zero-point energy contri¬ 
butions. Results are compared to experimental data extrapolated to 0 K. 
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T (K) Ci 2 (GPa) C' (GPa) Cm/C' 
0 (no ZPE) 151.4±0.2 72.7±0.3 1.44 
0 (ZPE) 148.01±0.5 70.0±0.4 1.46 
0 (Expt.)P 135.7 51.9 2.32 


TABLE IV. Ci 2 and C’ elastic constants, and Cm/C anisotropy ratio, derived from Tab. Ill with 
and without zero-point energy. Results are compared to experimental data extrapolated at 0 K. 


Errors are obtained according to propagation of uncertainties. 


C. Discussion 


The temperature dependence of the bulk modulus and of the elastic constants display an 
overall good agreement with the available experimental data, showing how lattice vibrations 
alone provide a robust description of the thermoelastic properties of the material, especially 
below the Debye temperature Qd- The agreement is still valid above Qd for the Cm-, 
that shows a near-linear behavior generally expected for metallic systems according to the 
semiempirical Varshni equation!^. On the other hand, approaching the Curie temperature 
(1043 K) from below, the results are not able to reproduce the anomalous non-linear softening 
which is observed in experimental B and Cn. 

According to previous work, e.g. based on a tight-binding approximation coupled to a 
single-site spin-fluctuation theory of band magnetismP^, effective spin-lattice couples mod¬ 
el^ as well as experimenters®, the origin of these anomalies is inherently related to magnetic 
fluctuations and, ultimately, to their influence on the free energy landscape (via modulation 
of the exchange couplings, conhgurational disorder and magnon-phonon interaction). In 
support of this conclusion, previous ab-initio papere^S® suggest that the electronic entropy 
and phonon-phonon a n h armonic effects beyond the quasiharmonic approximation play a 
minor role in determining the thermodynamics of the system below Tc. Given the strong 
indications of the pivotal role of magnetism in the description of thermoelastic properties 
of a-iron close to the Curie temperature, further ab-initio calculations taking into account 
magnetic disorder would be of paramount interest (see for instance Ref. I24|) . 

Focusing instead on the 0 K structural and elastic properties, we now discuss the possible 
origin of the mismatch between the calculated and experimental values. As we showed 
earlier, our calculated points are numerically accurate, and the errors associated with the 
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FIG. 12. Top panel: calculated adiabatic Cn elastic constant (blue solid line). Bottom panel: 
calculated adiabatic C 44 (blue solid line). Two sets of experimental data are reported on each plot 
- Expt.l (green circles) from Ref. 1501 and Expt.2 (yellow triangles) from Ref. [26l The calculated 
interval of confidence is displayed as a shaded area. As a guide to the eye, we also plot the elastic 
constants rigidly shifted (dotted line) and scaled (dashed line) to match the experimental values 
at zero temperature. 


fit are fairly small. As a consequence, the origin of the discrepancy can be ascribed (i) to 
the presence of magnetic domain walls, (ii) to the pseudopotential approximation, (iii) to 
the approximate XC functional. 

First, we inspected the possible effects on the equilibrium volume at 0 K due to the pres¬ 
ence of magnetic domain walls. We focus our attention to the collinear domain wall case, 
thus mimicking a magnetic distribution of domains in iron as two 8-atoms-thick ferromag- 
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FIG. 13. Thermal behavior of the C elastic constant calculated as a linear combination of B[T) and 
Cii(r) (blue continuous line). Two sets of experimental data are reported - Expt.l (green circles) 
from Ref. m and Expt.2 (yellow triangles) from Ref. EH The calculated interval of conhdence is 
displayed as a shaded area. As a guide to the eye, we also plot the elastic constants rigidly shifted 
(dotted line) and scaled (dashed line) to match the experimental values at zero temperature. 


netic strips with antiparallel magnetic moments repeated in the z direction throngh periodic 
bonndary conditions. The effect due the interfaces is to increase the lattice parameter of 
about 0.7% and to decrease the bulk modulus of about 9%. Since the density of domain 
walls in real materials is expected to be an order of magnitude lower, the effect on the 
lattice parameter should be rescaled accordingly, thus suggesting that domain walls affects 
only marginally the value of the equilibrium lattice parameter at 0 K. 


Next, we have observed that, for a given XC functional, details of the pseudopotential can 
have a large impact on the calculated quantities. For instance, in the case of PBE, the bulk 
modulus at 0 K ranges from 165 to 201 GPa if we consider pseudopotentials generated by 
different authors: ultrasoft or PAW, with 8 or 16 electrons or, even using the same electronic 
configuration but different version of pslibrarj^ (see Figs. Bin- 


As discussed in Sec. |III[ the pseudopotential used in this work was chosen as being closest 
in its equation of state and its magnetization as a function of volume to all-electron FLAPW 
calculation^^SEU result, the discrepancy at 0 K with respect to experiments found in 

this work and in all-electron calculations seems ascribable mainly to the exchange-correlation 
functional used. For this reason, we explored the effect of the XC functional on the 0 K 
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properties, keeping the pseudopotential generation scheme and parameters unchanged. We 
performed test calculations with the and PBEsoP^ functionals, and found that the 

disagreement with the experimental data is increased (see Figs. Bill- 

Eventually, we observe that the QHA thermal contribution to the energy landscape is 


almost linear (see Fig. 10) and does not contribute too much to the total curvature in 
the energy landscape. Moreover, its change in second derivative along with temperature is 
even smaller, and only marginally contributes to the temperature dependence of the total 
curvature of the energy landscape (its main effect is to shift the minimum of the free energy as 
a function of temperature). Therefore, we conclude hrst that the mismatch with experiments 
at hnite temperature is dominated by the 0 K static contribution discussed above. Second, 
that the temperature dependence of the elastic constants is driven, in hrst approximation, 
by the curvature of the 0 K energy landscape at the equilibrium expanded volumes. Our 
Ending suggests that one could try and employ more computationally expensive methods 
(such as BET+U+ J— hybrid functionald^, or to explore possible 

improvements in the description of the 0 K mechanical properties of a-iron, while thermal 
properties can be determined using lattice dynamics calculations performed with standard 
semi-local GGA functionals. 


V. CONCLUSIONS 

We have calculated the isothermal and adiabatic elastic constants of a-iron as a function 
of temperature from hrst-principles, using pseudopotential total energy calculations based 
on BET and lattice dynamics calculations based on BFPT, out of which we calculate free 
energies in the quasiharmonic approximation and hnite-temperature elastic constants from 
small strain deformations. Great care has been put in the verihcation of the pseudopo¬ 
tentials, and the validation of the results against experiments. Gommon semi-local BET 
functionals such as PBE reproduce only fairly elastic constants at zero temperature; on the 
other hand, their thermal behavior, originating from the changes in phonon dispersions upon 
crystal expansion, is very well described by the same functionals and in the quasiharmonic 
approximation, with a softening of the elastic constants and bulk modulus that is in ex¬ 
cellent agreement with experiments up to ©d and above. This work was supported by a 
grant from the Swiss National Supercomputing Gentre (GSGS) under project IB ch3. We 
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